# Loading packages
#library(tidyverse)
#install.packages("sf")
#install.packages("terra")
#install.packages("tmap")
library(tidyr)
library(sf) # simple features
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library (terra) # for raster
## terra 1.7.71
##
## Attaching package: 'terra'
## The following object is masked from 'package:tidyr':
##
## extract
library(tmap) # Thematic maps are geographical maps in which spatial data distributions are visualized
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(readr)
dat <- read_csv("../data/copepods_raw.csv")
## Rows: 5313 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): sample_time_utc, project, route, vessel, region
## dbl (6): silk_id, segment_no, latitude, longitude, meanlong, richness_raw
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Creating a figure that plots our spatial data, but isnt neccessarily spatially accurate
library(ggplot2)
ggplot(dat) +
aes(x = longitude, y = latitude, color = richness_raw) +
geom_point()

#So, now let’s look at the richness data
ggplot(dat, aes(x = latitude, y = richness_raw)) +
stat_smooth() +
geom_point()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

#turn our point data into a spatially referenced data frame
sdat <- st_as_sf(dat, coords = c("longitude", "latitude"),
crs = 4326)
#Coordinate reference systems
#Our copepod coordinates are long-lat, so we chose a common ‘one-size-fits-all’ GCS called WGS84 to define the crs using the EPSG code 4326.
crs4326 <- st_crs(4326)
crs4326 # look at the whole CRS
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
crs4326$Name # pull out just the name of the crs
## [1] "WGS 84"
crs4326$wkt # crs in well-known text format
## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
#Feature collection (points)
sdat
## Simple feature collection with 5313 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 89.6107 ymin: -65.2428 xmax: 174.335 ymax: -16.80253
## Geodetic CRS: WGS 84
## # A tibble: 5,313 × 10
## silk_id segment_no sample_time_utc project route vessel meanlong region
## * <dbl> <dbl> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 1 1 26/06/2009 22:08 AusCPR BRSY ANL Windar… 153. East
## 2 1 5 26/06/2009 23:12 AusCPR BRSY ANL Windar… 153. East
## 3 1 9 27/06/2009 0:17 AusCPR BRSY ANL Windar… 153. East
## 4 1 13 27/06/2009 1:22 AusCPR BRSY ANL Windar… 153. East
## 5 1 17 27/06/2009 2:26 AusCPR BRSY ANL Windar… 153. East
## 6 1 18 27/06/2009 2:43 AusCPR BRSY ANL Windar… 153. East
## 7 1 26 27/06/2009 4:52 AusCPR BRSY ANL Windar… 153. East
## 8 1 30 27/06/2009 5:57 AusCPR BRSY ANL Windar… 153. East
## 9 1 33 27/06/2009 6:45 AusCPR BRSY ANL Windar… 153. East
## 10 1 37 27/06/2009 7:50 AusCPR BRSY ANL Windar… 153. East
## # ℹ 5,303 more rows
## # ℹ 2 more variables: richness_raw <dbl>, geometry <POINT [°]>
# Cartography
# Plotting Richness
plot(sdat["richness_raw"])

plot(sdat)

#Thematic maps for communication
#Build and add on layers using using tmap
tm1<- tm_shape(sdat) +
tm_dots(col = "richness_raw")
# Saving our figure
tmap_save(tm1, filename = "Richness-map.png",
width = 600, height = 600)
## Map saved to /Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/code/Richness-map.png
## Resolution: 600 by 600 pixels
## Size: 2 by 2 inches (300 dpi)
# Mapping spatial polygons as layers
# First, load spatial data
aus <- st_read("../data/spatial-data/Aussie/Aussie.shp")
## Reading layer `Aussie' from data source
## `/Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/data/spatial-data/Aussie/Aussie.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 112.9211 ymin: -43.63192 xmax: 153.6389 ymax: -9.229614
## Geodetic CRS: WGS 84
shelf <- st_read("../data/spatial-data/aus_shelf/aus_shelf.shp")
## Reading layer `aus_shelf' from data source
## `/Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/data/spatial-data/aus_shelf/aus_shelf.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 112.2242 ymin: -44.1284 xmax: 153.8942 ymax: -8.8798
## Geodetic CRS: GRS 1980(IUGG, 1980)
# Mapping your polygons
tm_shape(shelf) +
tm_polygons()

# Layering copepod data points
tm_shape(shelf, bbox = sdat) +
tm_polygons() +
tm_shape(aus) +
tm_polygons() +
tm_shape(sdat) +
tm_dots()

#Exploring t_map
vignette('tmap-getstarted')
## starting httpd help server ... done
tm_shape(shelf) +
tm_polygons()

# Adding additional element, such as colors
tm_shape(shelf, bbox = sdat) +
tm_polygons() +
tm_shape(aus) +
tm_polygons() +
tm_shape(sdat) +
tm_dots(size = 0.01, col = "blue")

# Vignette Exercises
# World map
data("World") +
tm_shape(World) +
tm_polygons("HPI")

# Multiple shapes and layers
data(World, metro, rivers, land)
tmap_mode("plot")
## tmap mode set to plotting
## tmap mode set to plotting
tm_shape(land) +
tm_raster("elevation", palette = terrain.colors(10)) +
tm_shape(World) +
tm_borders("white", lwd = .5) +
tm_text("iso_a3", size = "AREA") +
tm_shape(metro) +
tm_symbols(col = "red", size = "pop2020", scale = .5) +
tm_legend(show = FALSE)
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

# Facets
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(World) +
tm_polygons(c("HPI", "economy")) +
tm_facets(sync = TRUE, ncol = 2)
#Splitting the spatial data
tmap_mode("plot")
## tmap mode set to plotting
## tmap mode set to plotting
data(NLD_muni)
tm1 <- tm_shape(NLD_muni) + tm_polygons("population", convert2density = TRUE)
tm2 <- tm_shape(NLD_muni) + tm_bubbles(size = "population")
tmap_arrange(tm1, tm2)
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.32. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

#Basemaps and overlay tile maps
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Legend for symbol sizes not available in view mode.